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ABSTRACT 


This  thesis  models  the  effects  in  a  fluid  medium  of  a  plane  wave  that  has  impinged  upon  a 
reinforced  plate.  The  wave  equation  for  pressure,  and  the  equation  of  a  thin  plate  combined  with 
other  equations  are  coupled  at  the  interface  between  the  fluid  and  the  thin  plate.  The  actual  modeling 
is  accomplished  in  a  fortran  computer  program  written  to  run  on  the  Naval  Postgraduate  School’s 
main  frame  computer.  The  program  uses  extensive  finite  differencing  on  a  domain,  assumed  to  be 
a  small  section  of  an  infinite  interface  between  the  fluid  and  the  plate,  to  simulate  the  deflection  of 
the  thin  plate  and  the  pressure  disturbances  in  the  fluid  medium.  To  accomplish  this,  each  of  the 
above  equations  will  be  scaled  or  nondimensionalized.  Additional  finite  differencing  is  explained 
which  covers  the  special  cases  for  the  side  boundaries  of  the  fluid  domain,  and  the  artificial  boundary 
created  to  model  infinity.  Different  beam  spacing  is  explored  for  its  effect  on  the  magnitude  of  the 
propagating  modes  of  the  scattered  pressure  wave. 
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I.  INTRODUCTION 


This  thesis  models  numerically,  the  effects  of  a  pressure  disturbance  incident  on 
a  fluid  loaded  reinforced  plate  using  finite  difference  methods.  The  program  produced 
explores  the  reaction  of  a  plane  wave  striking  the  plate  (hull  of  a  ship)  at  different 
frequencies,  for  different  plate  lengths  and  for  various  spacings  of  the  reinforcing  beams 
(or  members). 

Figure  1  shows  the  problem  being  modeled.  The  figure  shows  a  cross  section  of 
a  reinforced  flat  plate  with  water  on  the  top  side  of  it  and  a  vacuum  below  it.  For  this 
problem,  it  is  assumed  that  the  plate  is  infinite  in  the  x  direction  along  the  horizontal 
axis.  Since  the  plate  structure  is  periodic  in  space,  the  model  will  only  look  at  one 
section  of  this  infinite  plate.  The  behavior  at  other  sections  is  assumed  to  be  the  same. 
We  will  use  three  fundamental  equations  in  the  model.  The  first  equation  is  the  wave 
equation,  the  second  is  the  linear-inviscid-force  equation,  and  the  third  is  the  equation 
of  motion  for  a  thin  plate.  The  coupling  of  these  equations  at  the  fluid-plate  interface 
makes  the  problem  very  difficult  to  solve  analytically. 

Chapter  II  derives  the  relationships  among  these  three  equations  as  they  are  used 
in  the  model.  Assumptions  used  in  the  program  are  spelled  out  in  this  chapter.  Chapter 
II  also  includes  the  steps  necessary  to  scale  each  equation  to  facilitate  its  use  in  the 
computer  program. 
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Chapter  III  explains  the  operation  of  the  FORTRAN  computer  program  used  to 
model  the  problem.  The  program  takes  as  an  input  the  incident  and  specularly-reflected 
sound  waves  and  determines  the  scattered  pressure  due  to  the  fluid-structure  coupling. 
To  make  best  use  of  the  Naval  Postgraduate  School’s  mainframe  computer  the  main 
program  is  divided  into  a  number  of  smaller  subroutines.  These  are  also  described  in 
Chaptdr  III. 

Chapter  IV  discusses  the  results  of  the  program  runs.  It  also  reviews  the 
weaknesses  in  the  model  and  model  limitations. 


FLUID 


STIFFENER 

:  VAcuun 


Figure  1.  Cross  Section  of  Hull 
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II.  THEORY 


A.  THE  WAVE  EQUATION 

It  is  known  that  sound  travels  through  seawater,  or  any  fluid  medium,  as  a  pressure 
wave.  This  is  possible  because  seawater  is  a  compressible  fluid  which  is  affected  by 
pressure  disturbances.  To  approximate  the  movement  of  sound  in  a  fluid  medium  we 
have  used  the  linearized,  lossless  wave  equation  for  pressure  as  derived  in  Kinsler,  Frey, 
Coppens  and  Sanders  [Ref.  l:p.  104], 
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(2.1) 


where  Cf  is  the  fluid  sound  speed  and  p’^  is  the  total  pressure. 

The  pressure  wave  can  be  split  up  into  various  parts,  each  of  which  is  a  solution 
of  the  linear  wave  equation.  In  this  case  we  want  to  be  able  to  solve  for  pressures  that 
would  be  present  in  the  fluid  which  modify  the  pressure  resulting  from  specular  reflection 
of  an  incident  plane  wave  from  an  acoustically  hard  surface.  This  known  specularly 
reflected  pressure  wave  is  combined  with  the  computed  scattered  pressure  wave  to  make 
up  the  total  reflected  pressure  wave. 

Symbolically,  we  have: 
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(2.2) 


where  is  the  total  pressure,  p'  is  the  incident  pressure,  p"  is  the  specularly  reflected 
pressure,  and  p®  is  the  scattered  pressure.  The  total  radiated  pressure  is  a  superposition 
of  the  reflected  iuid  scattered  pressures. 

To  eliminate  the  specularly  reflected  wave  from  the  equations  of  motion  we  apply 
the  boundary  condition 

=0 

^  ay  ',.0 

In  this  case  the  normal  derivatives  of  the  incident  and  specularly  reflected 
waves  must  cancel  at  the  surface  of  the  plate.  The  vertical  coordinate  Y  is  equal  to  zero 
at  the  surface. 

We  assume  that  the  source  of  the  incident  pressure  is  far  from  the  surface  so  that 
a  plane  wave  approximation  is  adequate  to  model  the  incident  pressure,  p'. 

The  time  harmonic  incident  and  specularly  reflected  plane  waves  are  of  the  form: 


p\x,y,t)=€ 


(2.4) 
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and 


till 

jo( - -0 

I  V  -%i  M 


P  (x,y,i)=e 


(2.5) 


where 


COSQj 

-sin9^ 


and 


COS6y 

sin6. 


(2.6) 


(2.7) 


are  the  unit  direction  vectors  of  the  two  plane  waves. 

Both  the  incident  wave,  Equation  2.4  and  the  reflected  wave,  Equation  2.5  must 
satisfy  Equation  2.1.  Therefore,  because  of  superposition,  p*  must  also  satisfy  the  wave 
equation, 


^p- 


1  d^p 
cj  dt^ 


(2.8) 


B.  THE  LINEAR  BVVISCID  FORCE  EQUATION 

The  linear  inviscid  force  equation,  taken  from  Kinsler,  Frey,  Coppens  and  Sanders 
[Ref.  l:p.  104],  is  basically  the  inviscid  Navier  Stokes  equation  in  which  nonlinear  terms 
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have  been  neglected.  It  is  valid  for  acoustic  disturbances  of  small  amplitude,  and  can  be 
written 


where  Pf  is  the  constant  equilibrium  density  of  the  fluid,  and  u  is  the  velocity  of  the 
fluid. 

Since  we  assume  no  cavitation  of  the  fluid  or  penetration  of  the  fluid  at  the  fluid- 
plate  interface,  the  component  of  the  fluid  velocity  normal  to  the  interface  and  the  plate 
transverse  velocity  are  equal  at  the  fluid-plate  interface  at  all  times  t.  Letting  w  represent 
the  transverse  displacement  of  the  plate  (y  direction),  we  have 


=  at  y=0.  (2.10) 

dy 

Through  use  of  Equations  2.2  and  2.3  the  total  pressure  p^  can  be  replaced  by  the 
scattered  pressure  p®: 


dp^  , 

Of - =  --^—  at  y  =  0. 

^dt^  dy 


(2.11) 


C.  EQUATION  OF  MOTION  FOR  A  THIN  PLATE 

To  properly  understand  what  is  happening  on  and  within  the  thin  plate,  a  reference 
system  must  be  assigned.  In  Figure  2  we  see  that  the  vertical  plane  will  be  designated 
the  y-axis  in  which  the  thin  plate  experiences  transverse  displacement  w.  The  horizontal 
plane  will  be  the  x-axis  which  has  longitudinal  displacement  u.  To  set  up  our  equations. 
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we  assume  that  the  centerline  (axis)  of  the  beam  stays  perpendicular  to  its  cross  section 
during  bending. 


/ 
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- 21 - 

Figure  2  Plate  Reference  System 


As  put  forth  by  Kinsler,  Frey,  Coppens  and  Sanders  [Ref.  l:p.  58],  Hooke’s  law 
states  that  if  the  strain  is  small  the  stress  is  proportional  to  it.  The  stress  on  the  plate  is 
described  by  the  equation 


=  -y 


E 

1-v*  ax*' 


(2.12) 


where: 

is  the  stress  on  the  plate,  E  is  the  Young’s  modulus,  and  v  is  the  Poisson’s  ratio. 
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Consider  a  small  section  of  the  plate  of  area  AxAz  with  moments  and  forces  applied 
as  shown  in  Figure  3: 


BLOW-UP  Of  section  a-M 


Figure  3.  Forces  and  Moments 

Note  that  there  are  no  forces  or  moments  applied  which  have  z  dependence.  The 
motion  therefore  is  assumed  to  be  two-dimensional.  Since  the  plate  is  not  rotating,  a 
force  balance  is  performed  on  the  moments  and  moment  arm,  neglecting  terms  of  order 
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Ax^,  giving  us 


dM  r 

where: 

M  =  bending  moment  on  the  element, 

q  =  load  on  the  plate  element  (uniform  for  a  normally  incident  wave), 
and  F  =  flexural  force  applied  to  the  element. 

Equating  forces  in  the  y-direction  applied  to  the  section  of  plate  yields 

p  -^^(/jAxAz)  =  — AxAz-^AxAz, 
dx 

where: 

S  =  cross  sectional  area  of  the  plate, 
p,  =  density  of  beam  section, 
h  =  plate  thickness. 

Dividing  out  the  AxAz  terms  yields 


dF 

^  a/2 


Using  the  definition  of  the  moment 


h 

2 

_h 

~2 

where  the  moment  of  inertia  is  found  from 


El 

1-V2  dx^' 


(2.13) 


(2.14) 


(2.15) 


(2.16) 
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(2.17) 


h 

2 


l=fy^dy  = 

_h 

~2 

we  obtain  the  equation  for  deflection  of  a  plate 


El 

- +  p  A - 

i-v^ax^  ^  dt^ 


-<i- 


(2.18) 


The  loading  of  the  thin  plate  is  due  to  pressures  exerted  by  a  fluid  in  contact  with 
the  thin  plate.  Rewriting  Equation  2.18  including  the  pressure  terms  yields; 


in  which 


_  c^w  ,  ci^w  j 
D - +p^ - =  -p  , 

dx^  ^  dt^ 


(2.19) 


D=.  (2.20) 

12(1 -v^) 

is  the  flexural  rigidity  and  h  is  the  plate  thickness,  Junger  and  Feit  [Ref.  2:p.  194]. 


D.  THE  SCALED  EQUATIONS 
1.  The  Scaling  Factors 

In  order  to  use  the  three  governing  Equations  2.1,  2.9  and  2.19  in  the  finite 
difference  code,  they  must  be  scaled  or  nondimensionalized.  The  steps  are  as  follows. 
The  independent  variables  x  and  t  are  scaled  through  use  of  the  fundamental  length  L 
and  the  wave  frequency  w,  respectively. 
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and  y  =  —, 
L 


(2.21) 


x  = 


T 


where  x  and  y  are  the  scaled  spatial  coordinates,  and  L  is  the  one  half  the  distance 
between  stiffeners  measured  in  the  x-direction. 

The  equation  for  the  scaled  time  is 


t  =0)/, 

where  w  is  the  frequency  of  excitation  of  the  incident  plane  wave. 
We  then  have 


(2.22) 


a  _  a  _  1  a 

a.r  a(ZJf)  L  dx  * 


(2.23) 


and 


a  a  _  a 

ar  r  ■‘"aT'  (2-24) 

a(— ) 

w 

Additionally,  the  dependent  variables  w  and  p  must  be  scaled  as  follows.  For 
the  displacement  of  the  plate  (in  the  y-direction)  we  put 


-  w 


(2.25) 


where  w  is  the  scaled  vertical  displacement  of  the  plate  and  W  is  the  displacement 
scaling  factor.  For  pressure,  put 
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(2.26) 


where  p  is  the  scaled  pressure  and  =  Lw^pfW  is  the  pressure  scaling  factor. 

The  wave  numbers  for  waves  propagating  in  the  fluid,  and  along  the  plate  are 
given  as  follows, 


and 


(2.27) 


=r _ 

^  D  ^  ’ 


respectively. 

Additional  parameters  resulting  from  the  scaling  analysis  are 


(2.28) 


(2.29) 


where  fi  is  the  squared  ratio  of  fluid  to  plate  wave  numbers  and  Zy^,  the  coincidence 
frequency  of  the  plate,  is  defined  by 


We  also  define  the  dimensionless  quantity 


(2.30) 
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(2.31) 


where  e 

2. 


e=-^, 

P^^g’ 

is  a  dimensionless  quality  which  is  a  measure  of  the  fluid  loading  on  the  plate. 
The  Wave  Equation 

The  wave  equation  governing  the  scattered  fluid  pressure  was  found  to  be: 


dt^ 


(2.32) 


Scaling  the  above  equation  yields: 


Cj  L 


(2.33) 


which  reduces  to  the  desired  non-dimensional  form, 


-  -  r2/2- 


(2.34) 


3.  The  Linear  Inviscid  Force  Equation 
The  linear  inviscid  force  equation 


P/ 


dt^  dy 


in  scaled  form  reduces  to 


(2.36) 


(2.37) 
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4.  Thin  Plate  Equation 


The  dimensional  equation  of  flexural  rigidity  is: 


_  d*w  .  T\ 

By  substituting  in  the  scaling  factors  it  becomes: 

or,  after  more  simplification,  it  becomes 


(2.38) 


(2.39) 


w--  =  - 


'^kk*  _  _ 


(2.40) 


This  is  the  scaled  form  of  the  thin  plate  equation  that  will  be  used  in  the 
computer  program.  Note  that  at  y  =  0,  p*^  equals  p'.  Hence,  p’^  has  been  replaced  with 

p  +  2p’. 


5.  Stiffeners 

In  the  present  analysis  the  stiffeners  will  be  treated  in  an  idealized  way.  The 
end  conditions  we  will  assume  are  imposed  at  the  point  where  the  stiffeners  are  located 
will  be  equivalent  to  those  of  a  clamped  plate  for  which  the  boundary  conditions  are 
w  =  0  at  X  =  ±1,  and  w,  =  0  x  =  ±1. 


E.  THE  FLUID  DOMAIN 

As  stated  before,  it  is  assumed  the  plate  on  which  the  sound  wave  impinges  is 
infinite.  That  assumption  makes  it  possible  to  establish  the  following  conditions.  The 
boundary  conditions  on  the  section  of  plate  we  are  looking  at  will  be  periodic.  That  is 
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to  say,  for  normal  incidence  of  p',  displacement  on  the  plate  at  x  =  -1  will  be  the  same 
as  at  X  =  1.  For  arbitrary  angles  of  incidence,  the  pressures  at  x  =  ±1  differ  only  by 
a  scaling  factor  determined  by  the  angle  of  incidence.  Figure  4  and  the  following 
equations  will  show  how  the  time  harmonic  solution  is  dealt  with  in  the  problem. 


RADIAtiMS  Moper 

iK  vK 
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-v-h 
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LEFT 
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RiaHT 
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C  periodic) 

• 
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(periodic) 
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Figure  4.  Side  View  of  Fluid  Domain 


To  understand  how  the  fluid  pressure  is  disturbed,  the  wave  equation  is  considered 
for  which  a  prescribed  pressure  is  applied  from  the  fluid  to  the  surface  of  the  plate.  This 
uncoupled  problem  can  be  solved  analytically  by  the  method  of  separation  of  variables. 
The  scaled  wave  equation  is 
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(2.41) 


P  +/*  =  a^P  . 

*  XX  yy  **  ■*  tt 


Applying  periodic  boundary  conditions  results  in 


p(-l,y)=p(l,y)  and  p,( -1,50  =/>,(!.>')• 


(2.42) 


The  forced  pressure  at  the  plate  surface  is  given  by  the  boundary  conditions 


P(x,0,t)=8ix)e~*‘  where  g(-l)=g(l)  and  5,(-l)=g,(l), 

where: 

g(x)  =  prescribed  pressure  applied  to  the  fluid  by  the  plate. 

Assuming  a  time  harmonic  solution  with  equivalent  frequency  dependence  to  that 
of  the  forcing  function,  the  scattered  pressure  can  be  expressed  as 

By  substitution  of  Equation  2.44  into  Equation  2.41  we  obtain  the  reduced 
Hehmholtz  equation  for  </)(x,y), 

(2-45) 

with  boundary  conditions: 

4>(-1.50  =  4>(1,5^)  and  (|),(-l,y)=4>,(l,50  (2  46) 

and 
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4)(jr,0)  = 


(2.47) 


Henceforth  we  will  be  dropping  the  bars  from  all  equations,  giving  us  the 
following, 


(2.48) 

(t)(-l.y)=(t)(l.y). 

(2.49) 

(2.50) 

4>x(-i.y)=4>x(iy). 

(2.51) 

To  solve  for  4>  the  standard  method  of  "Separation  of  Variables"  is  applied  to  the 
equation.  First  we  write  as  a  product  of  functions  of  one  variable  each, 

<t)=X(;r)T(y).  (2-52) 

Upon  substitution  into  the  Helmholtz  equation  and  separation  of  variables,  we  find  a 
solution 

X  (x)-C  (2.53) 

'  II 

When  solved  for  Y,  two  answers  result: 

when  o>n7i,  (2-54) 


and 
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when  a  <  n  Tt  . 


(2.55) 


These  solutions  are  chosen  to  avoid  the  physically  implausible  solution  that  the  scattered 
pressure  as  y  approaches  infinity  will  blow  up  represent  an  incoming  wave. 

Combining  our  solutions  we  have 

••  u  _  «  _ 

<|)=5^<t:ii(x,y)='^  (2-56) 

n-O  #1=0  #i=A^+l 

Note  that  as  y  approaches  infinity  the  decaying  terms  become  exponentially  small.  Now 
from  our  boundary  condition  at  y  =  0,  q!„  can  be  found  such  that 

**  1  * 

g(x)=Y^  with  a^=— fg(x)e'"'^^dx.  (2.57) 

11=0  2^, 

It  is  important  to  note  that  in  the  numerical  solution  of  the  fully  coupled  problem, 
the  radiation  waves  must  be  properly  modeled  as  y  approaches  infinity.  The  numerical 
domain  cannot  be  semi-infinite  in  extent  in  the  y  direction.  It  must  be  truncated  at  some 
value  of  y,  which  we  call  yu,  where  a  boundary  condition  must  be  imposed  to  ensure 
proper  modeling  of  the  radiated  pressure. 


F.  THE  PERIODICITY  PROBLEM 

To  ensure  that  a  wave  with  any  angle  of  incidence  can  be  correctly  modeled  by  the 
computer  program  the  periodicity  must  be  worked  out.  If  we  were  to  have  a  wave 
normally  incident  on  the  plate  we  would  expect  a  solution  for  the  pressure  which  is 
periodic  in  x  to  be 
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(2.58) 

n-0 

To  account  for  an  arbitrary  angle  of  incidence,  Equation  2.58  needs  to  be  modified 
in  the  following  manner.  A  correction  factor  which  takes  into  account  the  phase  change 
at  X  =  -1  vice  x  =  1  must  be  added.  The  modified  equation  for  the  scattered  pressure 
is, 


n=Q 

where  6i  is  the  angle  of  propagation  of  the  incident  pressure  wave. 
Continuing  with  this  argument  one  can  write 

y^  =  -kcosBj-*-mi. 

Substitution  of  this  x  dependence  into  the  Helmholtz  equation  yields 

dy^ 


Substituting  these  relationships  back  into  our  equation  for  (f>  gives  us 

n--M  -«  N*\ 

where 


(2.59) 


(2.60) 


(2.61) 


(2.62) 
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and  in  which  the  flrst  sum  includes  all  values  of  n  for  which  k  >  7,  and  the  latter  two 
sums  account  for  the  remaining  modes. 


G.  PTNITE-DIFFERENCE  MODELING 
1.  The  Taylor  Series 

Refer  to  Numerical  Solution  of  Partial  Differential  Equations:  Finite 
Difference  Methods.  Chapter  I,  by  G.  D.  Smith  [Ref.  3]  for  a  more  complete  explanation 
of  finite  difference  methods. 

All  three  governing  equations  will  be  center  differenced  with  truncation  errors 
of  order  (Ax^,  At^). 

Finite  difference  approximations  for  u„  u„,  u^^^^  are 


(«x)r 


2h 


Oih\ 


(2.64) 


(2.65) 


and 


(«. 


‘  h* 


(2.65) 


where  Uj  =  u(Xi). 


2.  Differencing  of  the  Wave  Equation 


We  have  already  shown  that  the  wave  equation  is  of  the  following  form 


and  that  we  have  the  following  conditions, 

=  iAx,  yj  =  ;Ay,  r„  =  nAf, 


(2.67) 


(2.68) 


where  i,  j,  and  n  are  integers.  This  allows  us  to  describe  the  pressure  as  follows, 


=  P(x.,yj,t^)  and  Ax  =  Ay  =  /i. 

So  by  substituting  into  Equation  2.66  we  arrive  at  the  final  form. 


(2.69) 


ifi+i  At 


*2?," -Py (2.70) 


3.  Finite  Differencing  of  the  Linear  Inviscid  Force  Equation 

We  have  already  shown  that  the  linear  inviscid  force  equation  is  of  the 
following  form 


dp 

-w„  = 

"  dy 


(2.71) 


This  equation  is  applied  at  y  =  0  or  when  j  =  0  in  the  numerical  domain. 
The  displacement  w  does  not  depend  upon  y,  so  we  introduce  W“  =  W(Xj,  t„ )  and  write 
the  differenced  form  of  the  equation  as: 
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'  2Ay  '  ' 

This  equation  is  solved  for  the  P" which  is  exterior  to  the  numerical  domain. 
It  can  then  be  substituted  into  the  differenced  form  of  the  wave  equation  to  solve  for 
P“'*^’i  o  >  the  pressure  on  the  plate  fluid  interface  at  time  level  n+1  solving  for  Pj,.j  we 
have 


P"  =  p" 


(2.73) 


Note  that  is  first  needed  to  solve  for  P”  .,  which  means  a  prescribed  order  must  be 
followed  at  the  fluid  plate  interface.  First  W"^‘j  is  found,  the  P\.i,  and  finally 


4.  Finite  Differencing  the  Plate  Equation 

The  last  of  the  major  equations  that  we  must  put  in  the  form  of  finite 
differences  is  the  thin  plate  equation.  The  scaled  thin  plate  equation  is: 


*2?') .  P-'"') 

Differencing  this  equation  results  in  an  expression  for  the  displacement  W  at 
time  level  n  + 1  in  terms  of  the  preceding  two  time  levels  and  the  incident  plane  scattered 
pressure  at  time  level  n. 
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(2.75) 


^  (<2-^  ^i-2)  ^  Kx) 


Kl* 


h* 


i:^(p,V2p  \xfi^fi*iw^-wt-' 


5.  Finite  Differencing  Side  Boundaries 

When  differencing  the  wave  equation  at  the  side  boundaries  for  the  pressure 
we  have  a  similar  problem  to  that  found  in  applying  the  wave  equation  at  the  fluid  plate 
interface.  The  problem  is  that  there  are  values  of  the  pressure  needed  which  are  exterior 
to  the  numerical  domain.  Periodicity  will  be  used  as  an  additional  condition  to  eliminate 
this  problem.  We  have  the  scattered  pressure  in  the  form: 

<t>n(y>f)e‘”ne  (2.76) 

What  we  are  doing  here  is  taking  the  pressure  at  the  node  i  =  M  + 1,  x  =  Ax 
and  replacing  the  pressure  i  =  -M  +  1  and  x  =  -1  +Ax.  From  Equation  2.76  we  obtain 

(2-77) 

and  we  also  have: 

n.u‘E  (2,79) 


Now  we  can  make  the  following  two  statements: 
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(2.80) 


'“'E 


P-M*lJ=^  -*Axcos6,^«:cose,j^  (2.81) 


Combining  the  previous  statements  we  get: 


pH  -likcosBfnti 

An  equivalent  argument  used  at  the  left  boundary  yields: 


(2.82) 


p"  _-2i*cose,pn 


(2.83) 


6.  Finite  Differencing  the  Radiated  Boundary. 

The  final  edge  of  the  numerical  domain  is  the  artificial  boundary  which  must 
allow  free  passage  of  radiated  pressure  disturbances.  It  must  be  far  enough  away  so  that 
the  evanescent  modes  can  be  neglected.  To  see  how  the  boundary  condition  is  derived 
we  start  with  the  scattered  pressure  p’  and  try  to  derive  an  analytical  expression  for  the 
pressures  normal  direction  at  the  boundary,  as  is  done  in  Scandrett  and  Kriegsmann 
[Ref.  4].  This  type  of  radiation  boundary  condition  is  often  referred  to  as  a  non-local 
radiation  boundary  condition. 

From  previously  we  have 


N 

P\x,y/)=e''^  evanesent  modes),  (2.84) 

n^-M 

where 
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Y„  =  -ikcosQj-t^nn, 


(2.85) 


P.  =  A"-Yn 

and 

Additionally  in  our  problem  we  have 

k=k^L  and  ^/=— » 

but  as  y  approaches  infinity  p®  is  equivalent  to  only  the  radiated  modes, 
time  dependence  into  the  a„  coefficient  such  that, 

/i*-W 

and  in  the  following  steps  use  orthogonality  to  find  a„(t)e*^"’'.  We  have 
/.  "'■'p  *(w)&  =  ^  «,(r)« «’•” /« 


u[-ilxos6,»nn  -(-iJlcos0/+»iJi)] 
^u[nn-jnit] 

We  can  therefore  write 

where  |  =  dummy  variable  of  integration. 


(2.86) 

(2.87) 

(2.88) 

We  include  the 

(2.89) 

(2.90) 

(2.91) 
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(2.92) 


Noting  the  similarity  to  zero  in  the  following  equation 


h‘-m  dy  at 


Which  can  be  rewritten  as 


(2.93) 


dy 


(2.94) 


Remember,  these  equations  pertain  to  the  artificial  boundary  which  is  modeling  infinity. 
For  this  problem  y  is  equal  to  ye  at  the  artificial  boundary. 

When  we  start  to  finite  difference  Equation  2.94,  the  artificial  boundary  and 
use  the  trapezoidal  rule  of  integration  we  obtain: 


n—M  I— IX 


in  which  \ye  take; 


.  -  l  =  ±IX 

S,=  2 

1  l*±IX. 


(2.96) 


Rewriting  Equation  2.95  we  have 
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(2.97) 


^  Ml  ^  Mt 

2h 


h  “ 

*— E 


«,[E 


or 


pin  pm  /X 

Myy-i 


2/j 


*1, 

l‘-IX 


(2.98) 


Where  the  square  matrix  A  is  defined  as 


N 


Aj'— 5|E  P.‘ 


fyA-tj) 


Combining  the  above  equations  gives  us; 


(2.99) 


IX 


^:;K4rnV,^EA/(«-O-0 

l—K 

From  the  differenced  form  of  the  wave  equation  applied  at  j 
boundary)  we  have: 


(2.100) 

lY  (the  artificial 


Combining  Equations  2.100  and  2.101  and  eliminating  P“  jy+i  we  obtain: 
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(2.102) 


ft";'  -  -ft”;  ^-^[ft,:.jr’ft,:.jrt2fti..i 


iJY  ^  iJY  2  2 


A^ 


2  « 


«♦! 


kVihc  ‘■^  "  "  kv  " 


Now  we  need  to  put  this  equation  into  a  form  which  makes  it  possible  for  us 
to  see  the  matrixes  needed  to  solve  the  problem.  First  we  write; 


and 


p"'ii  *—y  A.pr^  ^ 

hVihx 


(2.103) 


We  define 


6  =/l 

\o  if  i*l 


(2.104) 


Aii  =  S.j+— A/ 


(2.105) 


B.,= 


(2.106) 


Equation  2.102  can  then  be  written  in  matrix  form  as 


=BP^'*  +  C 


(2.107) 


where 


28 


(2.108) 


Then  the  equation  for  the  pressure  vector  at  the  next  time  level  is  given  by: 


(2- 109) 

The  pressure  at  P"'^‘i.iY  is  then  assigned  by  taking  the  i*  component  of  the 
vector  P“^‘iY-  In  solving  the  matrix  Equation  2.107  we  use  the  UNPACK  library 
Gaussian  elimination  programs  (CGESO  and  CGESL). 
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III.  DESCRIPTION  OF  PROGRAM 


A.  GENERAL 

The  computer  program  written  for  this  thesis  simulates  the  effects  of  a  pressure 
disturbance  on  a  fluid  loaded  reinforced  thin  plate.  A  detailed  description  of  the  theory 
behind  the  program  can  be  found  in  Chapter  II.  A  pressure  wave  in  a  fluid  medium  is 
moving  towards  a  reinforced  plate  which  is  flat,  and  is  supported  by  stiffeners  which  are 
at  a  constant  spacing  on  the  non-fluid  side  of  the  plate.  The  non-fluid  side  is  assumed 
to  be  a  vacuum.  When  the  pressure  sound  wave  hits  the  plate  there  will  be  some  (a  very 
small  amount)  distortion  of  the  plate  in  the  vertical  direction.  This  displacement  will 
cause  a  disturbance  in  the  fluid  medium.  The  scalded  pressure  amplitude  for  different 
frequencies  and  beam  spacing  will  then  be  plotted  and  the  magnitude  of  the  propagating 
modes  will  be  recorded. 

The  program  stores  the  time  dependent  wave  equation  with  time  harmonic  forcing 
by  allowing  transient  effects  to  radiate  away  from  the  numerical  domain.  When  a  steady; 
state  has  ben  achieved  the  code  automatically  stops  the  calculations. 

B.  SUBROUTINES 

To  make  effective  use  of  the  Postgraduate  School’s  mainframe  computer,  the 
program  has  been  broken  down  into  a  number  of  smaller  programs  (i.e.  subroutines) 
which  can  be  more  efficiently  programmed  and  run. 
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The  main  program  is  used  to  coordinate  the  subroutines  and  to  pass  information 
from  one  subroutine  to  another.  The  following  are  descriptions  of  the  subroutines. 

1.  Subroutine  INIT 

The  INIT  subroutine  initializes  the  computer  program.  In  this  subroutine  we 
set  up  our  two  main  matrices.  The  first  of  these  is  the  "P"  matrix,  which  represents  the 
fluid  pressure  in  the  fluid  domain.  The  second  is  the  "W"  matrix,  this  represents  the 
displacement  of  the  plate  surface  which  bounds  the  fluid  medium. 

Although  we  have  the  ability  to  change  any  number  of  parameters  in  the 
program,  we  will  concentrate  on  two.  These  are  the  beam  spacing,  and  the  frequency 
of  the  pressure  wave  impinging  on  the  plate.  Both  of  these  values  are  assigned  in  this 
subroutine.  The  INIT  subroutine  sets  all  of  our  counters  and  subscripts  to  zero  at  this 
time  and  initializes  the  matrices  which  will  be  needed  to  model  the  artificial  boundary. 

2.  Subroutine  DOMAIN 

The  DOMAIN  subroutine  calculates  the  pressure  at  all  nodes  in  the  interior 
of  the  fluid  medium  at  a  subsequent  time  level  in  terms  of  the  preceding  pressure  values. 
The  subroutine  calculates  the  pressure  at  the  next  time  step  by  use  of  the  differenced 
wave  Equation  2.70. 

The  equation  above  must  be  modified  to  model  what  is  happening  at  the  side 
boundaries  of  the  fluid  domain.  As  was  stated  earlier,  our  model  is  looking  at  a  small 
section  of  an  "infinite"  fluid  plate  interface.  This  allows  us  to  use  the  pressure  in  the 
next  section  to  model  what  is  happening  in  the  section  in  which  we  are  interested.  This 
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is  necessary  since  otherwise  we  would  have  to  use  "undefined"  points  in  our  differenced 
wave  equation  at  the  edges  of  the  fluid  domain. 

3.  Subroutine  FLPL 

The  FLPL  subroutine  models  the  fluid-plate  interface  along  the  top  of  the 
reinforced  thin  plate.  Again,  with  finite  differencing  and  the  assumption  that  we  are  only 
looking  at  a  small  section  of  an  infinite  interface,  we  are  able  to  model  »he  effects  of  the 
pressure  wave  along  the  entire  fluid-plate  interface.  This  subroutine  is  also  where  we 
introduce  our  driving  function  into  the  computer  program.  This  function  is  what  starts 
things  moving  towards  a  steady  state.  Although  this  is  not  the  most  complicated 
subroutine  in  the  entire  program  it  is  the  one  which  has  proved  to  be  the  hardest  to  make 
work  as  predicted. 

4.  Subroutine  ARTF 

The  ARTF  subroutine  is  the  most  complicated  of  any  subroutine  in  the 
computer  program.  It  models  the  artificial  boundary  of  the  fluid  medium  at  infinity. 
Again,  this  is  accomplished  with  finite  differencing  and  matrix  operations. 

5.  Subroutine  ALPHA 

The  ALPHA  subroutine  calculates  the  magnitude  of  the  radiating  modes  of  the 
pressure  wave.  This  data  is  gathered  in  the  results  section  of  this  study. 

6.  Subroutines  SHIF  and  FINN 

The  SHIF  subroutine  moves  the  entire  problem  to  the  next  time  level.  This 
is  done  after  the  other  subroutines  have  completed  all  the  required  calculations  needed 
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to  satisfy  our  domain  and  all  of  our  boundary  conditions.  After  exiting  the  SHIP 
subroutine  the  entire  program  is  run  again  with  the  exception  of  the  INIT  subroutine. 

The  FINN  subroutine  is  a  check  to  make  sure  that  we  have  reached  a  steady  state 
situation  before  we  take  any  measurements  of  the  plate  displacement  or  of  any  fluid 
disturbances.  Figure  5  is  a  graphical  representation  of  the  displacement  convergence 
factor  approaching  zero.  This  was  arrived  at  with  a  assigned  frequency  of  3500  Hz  and 
L  =  1.0. 

7.  Subroutine  CLOSE 

The  CLOSE  subroutine  is  used  to  send  the  data  generated  in  other  parts  of  the 
main  program  to  the  plotting  devices. 

C.  PROGRAM  VERIFICATION 

Verification  of  the  program  was  accomplished  by  the  following  method.  We  solved 
the  waveguide  problem  (section  2E)  for  a  known  value  of  g(x).  We  let  g(x)  =  sin  ttx  and 
e'“.  This  was  accomplished  by  replacing  the  FLPT  subroutine  with  a  subroutine  named 
TEST.  As  before  the  program  was  run  until  a  steady  state  solution  was  reached.  These 
results  were  checked  against  analytical  solutions  and  found  to  compare  very  well.  For 
g(x)  =  sin  TX  we  graphically  reproduced  the  expected  results.  For  g(x)  =  e'”  we 
calculated  the  coefficients  a„  and  found  a,  -  1  and  all  other  a„’s  -  0.  The  small  error 
found  in  the  calculation  can  be  attributed  to  truncation  errors  in  the  finite  difference 
approximations.  These  results  match  with  analytical  calculations. 
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CONVERGENCE  VALUE 
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IV  RESULTS  AND  CONCLUSIONS 


A.  INITIAL  CONDITIONS 


To  properly  apply  the  computer  program  generated,  inputted  data,  or  conditions, 
must  be  within  the  tolerances  of  the  main  equations  used.  This  is  demonstrated  by  the 
use  of  .01  meters  as  the  thickness  of  the  thin  plate  and  a  frequency  range  from  500  Hz 
to  4000  Hz.  Other  inputted  data,  including  the  densities  of  steel  and  water  were  taken 
directly  from  Kinsler,  Frey,  Coppens  and  Sanders  [Ref.  l:pp.  461-463]. 

The  choice  of  the  Ax  and  Ay  terms  is  based  on  the  need  to  have  a  significant 
number  of  data  points  in  the  waveguide  domain  that  will  give  an  accurate  picture  of  what 
is  happening.  This  is  also  true  for  the  time  step  chosen,  we  need  to  make  sure  that  the 
At  was  small  enough  to  keep  the  problem  from  blowing  up.  A  Courant-Lewy-Friedrichs 
stability  condition  must  be  satisfied  for  the  domain  which  limits  the  size  of  At. 
Additionally  as  was  mentioned  before,  the  computer  program  was  run  for  the  same 
amount  of  time  for  each  trial,  this  time  period  is  long  enough  to  ensure  that  a  steady  state 
solution  was  being  achieved. 
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B.  PROGRAM  OUTPUT  DATA 


The  data  obtained  from  the  computer  program  demonstrates  the  following  major 
points.  When  the  frequency  is  increased  you  get  more  propagating  modes  from  each 
particular  plate  length.  When  L  is  changed  by  a  specific  amount  you  do  see  a  difference 
in  the  magnitude  and  shape  of  the  deformation  graphs. 

Our  first  table  shows  the  calculated  magnitudes  of  the  propagating  modes  of  the 
pressure  wave  at  an  assigned  frequency  of  750  Hz  at  different  beam  spacings. 


Table  I.  MAGNITUDE  OF  RADIATING  MODES. 


For  a  frequency  of  750  Hz 

with  L  =  0.8 

cvo  =  .467808127 

with  L  =  1.0 

ao  =  .484098911 

a,  =  .307484090 

with  L  =  1.2 

ao  =  .232444227 

a,  =  .367214561 

In  addition  to  the  data  in  Table  I,  we  have  in  Figures  6  through  8  a  graphical 
representation  of  the  scalded  pressure  amplitude.  The  pressure  was  measured  at  the 
artificial  boundary  ,  yn  at  a  frequency  of  750  Hz.  The  horizontal  axis  is  of  a  length  of 
2L  for  each  of  the  graphs  starting  at  x  =  -1  and  going  to  x  =  1 .  The  numbering  on  the 
horizontal  axis  represents  each  of  the  41  pressure  nodes  used  in  the  finite  difference 
code. 
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Figure  6.  Sample  Pressure  Graph 
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SCALED  PRESSURE  AMPURJDE 
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SCALED  PRESSURE  AMPLITUDE 
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Table  II  is  a  list  of  the  magnitudes  of  the  q!„’s  for  program  runs  with  a  frequency 


of  3500  Hz. 


Table  H.  MAGNITUDE  OF  RADIATING  MODES. 


For  a  frequency  of  3500  Hz 

with  L  =  0.8  ao  =  .451181054 

a,  =  .323620617 

a,  =  .088650584 

with  L  =  1.0  tto  =  .365310729 

a,  =  .260265648 

a,  =  .062102031 

a3  =  .060141488 

a4  =  .076183497 

with  L  =  1.2  ao  =  .252027571 

a,  =  .165709317 

a2  =  .047427550 

a3  =  .046469088 

a4  =  .049672558 

as  =  .043821129 

a^  =  .067380726 


From  the  data  above  we  observe  that  at  a  frequency  of  3500  Hz  there  is  a 
substantial  change  in  the  magnitude  of  the  propagating  modes  as  L  is  increased  from  0.8 
to  1.2.  It  can  be  seen  from  other  runs  that  the  higher  the  frequency  the  larger  the 
differences  in  the  magnitudes  of  the  propagating  modes. 

Figures  8  through  10  show  the  scaled  pressure  amplitude  data  obtained  from  the 
program  runs  with  an  assigned  frequency  of  3500  Hz.  It  is  observed  that  changing  of 
the  length  L  on  which  the  Ax  is  scaled  has  an  effect  on  not  only  the  amplitude  of  the 
pressure  wave  but  also  their  shape. 
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L  =  0.8  FREQUENCY  =  3500  HZ 


HORIZONTAL  POSITION 


LEGEND 

□  =  FLUID  PRESSURE 


Figure  9.  Sample  Pressure  Graph 


PRESSURE  AMPLITUDE 
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SCALED  PRESSURE  AMPLITUDE 


==  1.2  FREQUENCY  =  3500  HZ 


LEGEND 

□  =  FLUID  PRESSURE 


Figure  11.  Sample  Pressure  Graph 
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Table  III  is  a  list  smaple  energy  calculations  taken  at  three  frequencies  inputted  into 
the  computer  program.  The  data  presented  demonstrates  that  for  theses  frequencies  and 
lengths  run  in  the  computer  program  that  the  fundamental  mode  contains  the  largest 
amount  of  energy.  Also  evident  is  that  the  energy  contained  in  the  fundamental  mode 
starts  becomes  more  constant  as  the  frequency  is  increased  for  different  lengths  of  L. 
This  data  is  in  agreement  with  the  theory  of  sound  propagation. 
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Table  HI  ENERGY  CALCULATIONS 


Percent  of  Enerpv  in 

Percent  of  Energv  in 

Pronasating  Mode 

Prooaeatine  Mode 

For  a  frequency  of  750  Hz 

For  a  frequency  of  3500  Hz 

with  L  =  0.8  Eo  =  100% 

with  L  =  0.8  Eo  =  66% 

E,  =  32% 

with  L  =  1.0  Eo  =  99% 

E^  =  2% 

E,  =  1% 

with  L  =  1.2  Eo  =  77% 

with  L  =  1.0  Eo  =  65% 

E,  =  23% 

E,  =  32% 

Ej  =  2% 

For  a  frequency  of  1750  Hz 

Ej  =  1% 

E4  =  * 

with  L  =  0.8  Eo  =  77% 

E,  =  23% 

with  L  =  1.2  Eo  =  64% 

E,  =  27% 

with  L  =  1.0  Eo  =  68% 

E2  =  2% 

E,  =  31% 

Ej  =  2% 

Ej  =  1% 

E4  =  2% 

Ej  =  1% 

with  L  =  1.2  Eo  =  64% 

E,  =  2% 

E,  =  32% 

Ej  =  3% 

Ej  =  1% 

*  denotes  much  less  then  1  % 
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Table  IV  shows  the  magnitude  of  the  primary  propagating  mode  for  each  of  the 


different  frequencies  and  lengths  inputted  into  the  computer  program. 


Table  IV  MAGNITUDES  OF  PRIMARY  MODES 


Frgqygnpy 

750  Hz 

Magnitude  of  ofo’s 

L  =  0.8  L  =  1.0 

.467808127  .484098911 

L  =  1.2 

.367214561 

1000  Hz 

.6182562721 

.444052100 

.262355387 

1250  Hz 

.553467214 

.457091272 

.343211770 

1500  Hz 

.498957574 

.415586054 

.320898831 

1750  Hz 

.404151559 

.431289613 

.338643909 

2000  Hz 

.486700475 

.433078885 

.316784620 

2250  Hz 

.394520938 

.400619507 

.252146780 

2500  Hz 

.385663450 

.392564356 

.267684340 

2750  Hz 

.385663450 

.392564354 

.267684340 

3000  Hz 

.420025572 

.353338659 

.234915197 

3250  Hz 

.369634628 

.226448834 

.194508076 

3500  Hz 

.451181054 

.365310729 

.252027571 

3750  Hz 

.446322470 

.330691218 

.226533830 
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The  data  obtained  from  the  computer  program  and  represented  in  Table  III  shows 
that  there  is  always  a  difference  in  the  magnitude  of  the  fundamental  propagating  mode 
for  changes  in  L  at  all  frequencies  comprised  within  the  computer  runs.  The  difference 
in  the  magnitude  of  the  primary  propagating  mode  between  the  runs  with  a  length  "L" 
of  1.2  and  the  runs  with  a  length  L  of  1  are  greater  then  the  difference  between  the  runs 
with  lengths  of  0.8  and  1.0.  In  the  frequency  range  covered  when  L  is  increased  from 
1.0  to  1.2  there  is  always  a  drop  in  the  magnitude  of  the  fundamental  propagating  mode. 
This  is  not  always  true  when  L  is  decreased  to  0.8  from  1.0.  Although  the  magnitude 
of  these  increases  are  smaller  then  the  magnitude  of  the  decreases  this  demonstrates  that 
by  increasing  the  spacing  of  the  reinforcing  beams  will  not  always  result  in  a  decrease 
in  the  amplitude  of  the  fundamental  propagating  mode  in  this  frequency  range.  From  this 
result  we  can  conclude  that  additional  factors  play  an  important  role  in  determining  the 
magnitude  of  the  propagating  modes. 

Because  the  model  is  allowed  to  run  until  a  steady  state  solution  is  reached  there 
is  no  appreciable  difference  in  the  absolute  pressure  measured  at  the  artificial  boundary 
for  any  of  the  frequencies  observed.  However,  it  is  evident  that  as  the  length  of  the  plate 
"L"  is  increased  there  is  a  slight  drop  in  the  pressure  measurement.  This  is  true  for  all 
frequencies.  Figure  12  is  a  sample  plate  deformation  graph  for  L  =  1.0  and  a  frequency 
of  1000  Hz. 

In  conclusion,  the  computer  program  produces  results  as  expected  for  the  different 
frequencies  of  the  plane  wave  and  lengths  of  the  thin  plate.  The  superpositioning  of  the 
reflected  plane  wave  has  showed  that  as  L  is  increased  there  is  a  larger  number  of 
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propagating  modes  accompanied  by  a  change  in  the  amplitude  of  each  mode.  Additional 
results  could  be  explored  for  different  angles  of  deflection  as  well  as  different  plate 
thicknesses.  This  will  require  modification  of  the  computer  program  as  it  stands  to 
gather  data  for  additional  angles  and  thicknesses. 
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SCALED  DISPLACEMENT  OF  PLATE 


SCALED  LENGTH  OF  1.0 


LEGEND 

D  =  DISP  AT  1000  HZ 


Figure  12.  Deformation  Graph 
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APPENDIX:  COMPUTER  PROGRAM 


FILEi  280CT  FORTRAN  A1 


CXXXXXKKXXXXXKXKKXKXKXKKXXXXXXKXKXXXXKKXXXKXafXXXXXXXXKKXXXXXX 

.  THIS  PROGRAM  SIMULATES  THE  DISTRUBANCE  IN  A  FLUID 
AFTER  A  SOUND  HAVE  HAS  HIT  A  PLATE. 

P  =  PRESSURE  (EXCESS)  D  =  FLEXURAL  RIGIDITY 

X  =  H0RI20NTIAL  POSITION  Y  =  VERTICAL  POSITION 

W  =  DEFORMATION  KF  =  FLUID  WAVE  NUMBER 

SPD  =  WAVE  SPEED  KP  =  PLATE  WAVE  NUMBER 

FREQ  =  FREQUENCY  OMEGA  =  RATIO  OF  (KF/KP)»»2 

COUNTERS!  I  J  N 

THE  VALUES  IN  THE  PARAMETER  STATEMENTS  ARE  CONSTANTS 
PARAMETER  (IX  =  20,  lY  =  80) 

PARAMETER  (H  =  .025,  DELTAT  =  .0325  ) 

PARAMETER  (NDIM  =  2XIX+1) 

THE  FOLLOWING  ARE  COMPLEX  ARRAYS 
C0MPLEXK16  P(  -IXs  IX,  0i  lY,  3),  W(  -IXi  IX,  3) 
COMPLEX  LB,  FACR 

COMPLEX  A(  -IXt  IX,  -IXt  IX) 

COMPLEX  AA(  NDIM,  NDIM) 

COMPLEX  BB(  NDIM,  NDIM) 

COMPLEX  C(  -IXt  IX) 

COMPLEX  Z(  NDIM  ) 

C0MPLEXX16  FI 

COMMON/SS/P,  W 

COMMON/TT/KP,KF,INFA,  THETAI,L,  X,  PI, 

I  EISP,  OMEGA,  T 

COMMON/CC/LB,  FACR 
COMMON/EE/WGTWO,  WGTHE,  WGTOT,  CON 
COMMON/PT/PGTWO,  PGTHE,  PGTOT 
COMMON/HH/MINN,  MAXX 
COMMON/II/AA,  BB 
COMMON/KK/IPVT 
COMMON/MM/NN 
COMMON/OO/C 
COMMON/ FF/ A 
COMMON/UP/AL 
COMMON/RN/RUNS 
COMMON/FEL/U,V,TP 
COMMON/SN/FREQ,  SPD 
COMMON/GF/WGGRF, PGGRF 

INTEGER  RUNS,  IPVT(NDIM),  SAM 
PARAMETER  (SAM  -  1000) 

REAL  KF,  KP,  INFA,  EISP,  OMEGA,  L,PI 

REALX8  WGTWO, WGTHE, WGTOT, CON 

REALX8  PGTWO, PGTHE, PGTOT 

REAL  WGGRF(  SAM  ),  PGGRF(  SAM  ) 

REAL  THETAI,  FREQ,  SPD 

REAL  X,  T 

PI  =  3.141592654 
CON  =  .015D0 
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WGTWO  =  0.0 DO 
WGTHE  =  O.ODO 
PGTWO  =  O.ODO 
PGTHE  =  O.ODO 


CALL  XUFLOM  (0) 
CALL  INIT 
RUNS  =  0 
20  CALL  FLPL 
CALL  DOMAIN 

RUNS  =  RUNS  +  1 


C 


C 


CALL  ARTF 

IF  (RUNS  .GT.  1999)  THEN 

GOTO  60 

ELSE 

GOTO  61 

ENDIF 

60  CALL  FINN 

61  T  =  T  +  DELTAT 
CALL  SHIF 

IF  (RUNS  .EQ.  3000)  THEN 

GOTO  80 

ELSE 

GOTO  20 

ENDIF 

80  CALL  ALPHA 
CALL  CLOSE 


WRITE(2,767)THETAI,KF 

WRITE(2)768)0MEGA,EISp''^^'*'^^''  '  *'^12. 6)  , 

A87  FORMAT(i)(,'  H  =  '/F12  6»3X#*  L  =  •  F6 
WRITE(2,735)FREQ;sPD  >•  ''^^.3) 

WRITE(2!629)pP  '  »P12.6,3X,  •SOUND  SPEED  =  •,F12.6) 

629  FORMAT(ix,'  PI  =  •,F12.7) 

WRITE(2,919) 

919  FORMATdX, ‘FACR  =  •) 

WRITE(2,*)  FACR 
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WR1TE(2  222) 

222  FORMAKIX, 'FINAL  DATA* ,//, IX,  •  X 

1  •  W  P  •) 

DO  333  I  =  -IX,  IX,  I 
XX  =  IXH 

333  MRITE(2,666)XX,CDABS(M(I,2)),CDABS(P(I,IY,3)) 

666  FORMAT ( I X,5FI 0.5) 

WRITE(2,223)RUNS 

223*  FORMATdX, 'RUNS  =  ',17) 

WRITE(2,678) 

678  FORMATdX, 'BELOW  IS  THE  FINAL  WGTOT,  AND  PGTOT') 
WRITE(2,X)WGT0T 
WRITE(2,*)PGTOT 
STOP 
END 


SUBROUTINE  INIT 


THIS  SETS  UP  THE  INITIAL  CONDITIONS  OF  THE 
PROBLEM  INCLUDING  THE  DOMAIN  AND  THE 
ARTIFICIAL  BOUNDARY 


C 


THE  VALUES  IN  THE  PARAMETER  STATEMENTS  ARE  CONSTANTS 
PARAMETER  dX  =  20,  lY  =  80) 

PARAMETER  (H  =  .025,  DELTAT  =  .0325) 

PARAMETER  (NDIM  =  2KIX+I) 

THE  FOLLOWING  ARE  COMPLEX  ARRAYS 
COMPLEXXI6  P<  -IXt  IX,  Ot  lY,  3),  WC  -IX.  IX,  3) 
COMPLEX  LB,  FACR 

COMPLEX  A<  -IX. IX,  -IX.  IX) 

COMPLEX  AA(  NDIM  ,  NDIM  ),  BBC  NDIM  ,  NDIM  ) 

COMPLEX  Z(  NDIM  ) 

COMMON/SS/P,  W 

COMMON/TT/KP,KF,INFA,  THETAI,L,  X,  PI, 


FILE.  280CT  FORTRAN  A1 


C 


1  EISP,  OMEGA,  T 

COMMON/CC/LB,  FACR 
COMMON/ FF/ A 

COMMON/ GG/GAMMA,  BETTA 
COMMON/HH/MINN,  MAXX 
COMMON/II/AA,  BB 
COMMON/KK/IPVT 
COMMON/FEL/U,V,TP 
COMMON/SN/FREQ,  SPD 


REAL 

REAL 

REAL 

REAL 

REAL 

INTEGER  I 

INTEGER 

INTEGER 


THETAI,  INFA,  KF,  EISP,  OMEGA, 
ART,  K,  XI,  XIL,  RCOND 
GAMMA(-A0«  60),  BETTA  (-60.60), 
SPD,  FREQ 
X,  T 

,  J,  N,  LJ,  NN,  U,  V,  TP 
NMODMIN,  NMODMAX,  MINN,  MAXX 
IPVTCNDIM) 


TIN, 

PI 


L,  KP 
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TP  =  41 

LB  =  CMPLX(0.0,  1.0) 

FACR  =  CMPLXCl.O,  0.0) 

SPD  =  1500. ODO 
FREQ  =  3750. ODO 
N  =  0 
1  =  0 
J  =  0 
NN  =  1 
L  =  1.2D0 

THETAI  =  PI/2. ODO 

T  =  DELTA! 

NMODMIN  =  0 
NMOOMAX  =  0 

KF  =  (2D0KPI»FREQ*L)/SPD 

,  GT  =  ((195000000000. 0*( .05»«3))/12(l-( .28XX2)) 
KP  =  ((((FREQ*2*PI)»*2)»7700*H)/FR)*x.25 
KP  =  40.7534975147D0 
OMEGA  =  (KF/KP)»X2 

INFA  =  (DELTATXX2)/((KFX«2)K(L»»2)X(H»*2)) 

DO  12  I  =  -IX,IX,1 
DO  11  J  =  0,IY,1 

P(I,J,1)  =  CMPLX(0.0,  0.0) 

P(I,J,2)  =  CMPLX(0.0,  0.0) 

P(1,J,3)  =  CMPLX(0.0,  0.0) 

11  CONTINUE 

12  CONTINUE 

DO  21  I  =  -IX,  IX,  1 

W(I,1)  =  CMPLX(0.0,  0.0) 

W(I,2)  =  CMPLX(0.0,  0.0) 

W(I,3)  =  CMPLX(0.0,  0.0) 

21  CONTINUE 

THIS  SETS  UP  THE  ARTIFICAL  BOUNDRAY 
K  =  kF^L 

TIN  =  -KXCOS(THETAI) 

DO  15  N  =  -40,  40,  1 

GAMMA(N)  =  TIN  +  PIKN 

IF  (  ABS(GAMMA(N))  .LE.  K)  THEN 
NMODMIN  =  MIN(N,  NMODMIN) 

NMODMAX  =  MAX(N,  NMODMAX) 

ENDIF 

15  CONTINUE 

MINN  =  NMODMIN 
MAXX  =  NMODMAX 


DO  18  N  =  MINN,  MAXX,  1 

BETTA(N)  =  SQRT((KX*2)  -  (GAMMA(N)*K2) ) 
18  CONTINUE 

DO  78  I  =  -IX,  IX,  I 
DO  79  LJ  =  -IX,  IX,  1 

A(I,LJ)  =  CMPLX(0, 0,0.0) 

79  CONTINUE 
78  CONTINUE 

ART  =  (2*(HK»2))/(4»DELTAT) 

DO  300  I  =  -IX,  IX,  1 
DO  301  LJ  =  -IX,  IX,  1 
XI  =  I*H 
yTi  s  I  ivu 

IF  (ABS(LJ)  .EQ.  IX)  THEM 
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DEL  =  .5 
ELSE 
DEL  =  1 
ENDIF 

DO  302  N  =  MINN,  MAXX,  1 
A(I,LJ)  =  A(I,LJ)  +  ARTXDEL 

1  XBETTACN)«CEXP(LBXGAMMA(N)»{XI-XIL)) 

302  CONTINUE 
301  CONTINUE 
300  CONTINUE 
C 

DO  420  U  =  1,  NDIM,  1 
>  DO  421  V  =  1,  NDIM,  1 
I  =  -IX+U-1 
LJ  =  -IX+V-1 

AA(U,V)  =  INFA)((A(I,LJ)) 

BB(U,V)  =  INFAK(A(I,LJ)) 

IF  (I  .EQ.  LJ)  THEN 

AA(U,V)  =  AA(U,V)  +  1.0 
BB(U,V)  =  BB(U,V)  -  1.0 

ENDIF 

421  CONTINUE 
420  CONTINUE 

CALL  CGECO  (AA, NDIM, NDIM,  IPVT,  RCOND,Z) 

RETURN 

END 

C 

SUBROUTINE  FLPL 

THIS  SUBROUTINE  DEFINES  THE  FLUID-PLATE  INTERFACE 
THIS  IS  WHERE  Y=0  ALONG  THE  TOP  OF  THE  PLATE. 

THE  VALUES  IN  THE  PARAMETER  STATEMENTS  ARE  CONSTANTS 
PARAMETER  (IX  =  20,  lY  =  80) 

PARAMETER  (H  =  .025,  DELTAT  =  .0325) 

PARAMETER  (NDIM  =  2*IX+1) 

C  THE  FOLLOWING  ARE  COMPLEX  ARRAYS 

C0MPLEXX16  P(  -IXs  IX,  0«  lY,  3),  W(  -IX:  IX,  3) 
COMPLEX  LB,  FACR,  PTEMP 

C0MPLEXX16  FI 
COMMON/SS^P-  M 

COMMON/TT/RP,KF,INrA,  THE1AI,L,  X,  PI, 

1  EISP,  OMEGA,  T 

COMMON/CC/LB,  FACR 
COMMON/RN/RUNS 
C 

REAL  KP,  L,  OMEGA,  EISP,  KF,  INFA,  PI,  THETAI 
REAL  X,  T 

REALK4  SQ,  Si/1,  BW1,BW2,BW3,BW4 
INTEGER  I,  J,  RUNS 


EISP  =  0.13101864566D0 
C 

J  =  0 

SQ  =  (DELTAT/(H»KF))X»2 
SW  =  (DELTAT/(H*KP)KX2)x*2 
BWl  =  2-6 XSW 
BW2  -  4XSW 

BW3  =  EISPXDELTATXK2KKF/0NEGA 
BW4  =  -SW 

A  =  FL0AT(IX)/FL0AT(2XIX) 
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THIS  LOOP  COMPUTES  THE  INTERIOR  NODES 

DO  10  I  =  -IX  +  2,  IX  -2,  1 
X  =  1*H 

W(I,3)  =  BW1»W(I,2)  -W(I,1)  +  BW2*(M(I+1,2)+W(I-1,2)) 

1  +BW4*(W<I+2,2)+W(I-2,2)) 

2  +BW3»(P(I,0,2)+  F1(X,T)) 

PTEMP  =  -(2*H/DELTAT»*2)«(W(I,3)-2»W(I,2)+W(I,1))  +P(I,1»2) 

P(I,J,3)  =  INFAx(P(I+l, J,2)+P(I-1,J,2)-4»P{I, J,2)+P(I, J+l,2)+ 

1  PTEMP  )+2»P(I,J,2)  -PCI,J,1) 

10  CONTINUE 

THESE  LOOPS  COMPUTE  THE  NODES  AT  -19  AND  19 

I  =  -IX+1 
X  =  (-IX+1)XH 

WCI,3)  =  BW1XW(I,2)-W(I,1)  +BW2K(WCI+1,2)+W(I-1,2)) 

1  +BW4x(W(I+2,2)+W(I,2)) 

2  +DW3X(P(I,0,2)+  F1(X,T)) 

PTEMP  =  -(2XH/DELTATXX2)«(H(I,3)-2»W(I,2)+W(I,1))  +P(I,1»2) 

P(I,J,3)  =  1NFA»(P(I+1,J,2)+P(I-1,J,2)-4XP(I,J,2)+P(I,J+1,2)+ 

1  PTEMP  )+2xP(I,J,2)  -P(I,J,1) 


I  =  IX  -1 
X  =  (IX-I)XH 

M(I,3)  =  BH1XW(I,2)  -  W(I,1)  +  BW2K(W(I+1,2)+W(I-1,2)) 

1  +BH4X(H(I,2)+W(I-2,2)) 

2  +BW3x<P(I,0,2)+  F1(X,T)) 

PTEMP  =  -(2XH/DELTAT»X2)K(M(I,3)-2»W(I,2)+W(I,1))  +P(I,1,2) 
P(I,J,3)  =  INFAX(P(I+1,J,2)+P(I-1,J,2)-4XP(I,J,2)+P(I,J+1,2)+ 

1  PTEMP  )+2xp(I,J,2)  -P(I,J,1) 

THESE  ARE  THE  END  POINTS 

FACR=CEXP<2XLBXLXKFXC0S(THETAI)) 

I  =  -IX 
J  =  0 

W(I,3)  =  CMPLX(0. 0,0.0) 

PTEMP  =  P<I,1,2) 

P(I,J,3)  =  INFAX(P(I+1,J,2)+FACRXP(-I-1,J,2)-4XP(I,J,2. 

1  +P(I,J+1,2)+PTEMP  )+2xP(I,J,2)  -P(I,J,1) 

FACR=CEXP(-2xLBXLxkFXC0S(THETAI) ) 

I  =  IX 
J  =  0 

W(I,3)  =  CMPLX(0. 0,0.0) 

PTEMP  =  PCI, 1,2) 

P(I,J,3)  =  INFAXCFACRXPC-I+1,J,2)+P(I-1,J,2)-4XPCI,J,2) 

1  +P(I, J+1,2)+PTEMP  )+2xPCI,J,2)  -P(I,J,1) 

RETURN 

END 
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FILE.  280CT 


FORTRAN  A1 


SUBROUTINE  DOMAIN 

THE  VALUES  IN  THE  PARAMETER  STATEMENTS  ARE  CONSTANTS 
PARAMETER  CIX  =  20,  lY  =  80) 

PARAMETER  (H  =  .025,  DELTAT  =  .0325) 

PARAMETER  (NDIM  =  2»IX+1) 

THE  FOLLOWING  ARE  COMPLEX  ARRAYS 
COMPLEXX16  P(  -IX.  IX,  0.  lY,  3),  W(  -IX.  IX,  3) 
COMPLEX  LB,  FACR 

COMMON/SS/P,  W 

COMMON.'TT/KP,  KF,  INFA,  THETAI,  L,  X,  PI, 

1  EISP,  OMEGA  ,T 

COMMON/CC/LB,  FACR 

REAL  KP,  RF,  INFA,  THETAI,  L,  PI,  EISP,  OMEGA 
REAL  X,  T 
INTEGER  I,J 

DO  22  I  =  1-IX,  IX-1,1 
DO  21  J  =  1,  IY-1,  1 

P(I,J,3)  =  INFA«(P(I+1,J,2)+P(I-1,J,2)- 

1  4XP(I, J,2)+P(I,J+l,2)+P(I,J-l,2))+ 

2  2XP(I, J,2)-P(I,J,1) 

21  CONTINUE 

22  CONTINUE 

LB  =  CMPLX  (0.0,  1,0) 

LA  =  LXKFXCOS(THETAI) 

FACR  =  CEXP<(2)X(LBxlA)) 

I  »  -IX 

DO  41  J  =  1,  IY-1  ,  1 

P(I,J,3)  =  INFAX(P(I+1,J,2)+  FACRXP(-I-1,J,2) 

1  -4XP(I, J,2)+P(I, J+1,2)+P(I, J-l,2))+ 

2  2XP(I, J,2)-P(I,J,1) 

41  CONTINUE 

LB  =  CMPLX  (0.0,  1.0) 

RA  =  LXKFXCOS(THETAI) 

FACR  =  CEXP((-2)X(LBXRA)) 

I  =  IX 

DO  31  J  =  1,  IY-1,1 

P(I,J,3)  =  INFAX(FACRXP(-I+1,J,2)+  P(I-1,J,2) 

1  -4XP(I,J,2)+P(I,J+1,2)+P(I,J-1,2))+ 

2  2XP(I,U,2)-P(I,J,1) 

31  CONTINUE 

RETURN 

END 


56 


oooo 


C 
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SUBROUTINE  ARTF 

THE  VALUES  IN  THE  PARAMETER  STATEMENTS  ARE  CONSTANTS 
PARAMETER  (IX  =  20,  lY  =  fO) 

PARAMETER  (H  =  .025,  DELTAT  =  .0325) 

PARAMETER  (NDIM  =  2XIX+1) 

THE  FOLLOWING  ARE  COMPLEX  ARRAYS 


C0MPLEXX16 
COMPLEX 
COMPLEX 
COMPLEX 
COMPLEX 
COMPLEX 
COMPLEX 

COMMON/SS/P,  W 
COMMON/TT/KP,KF,INFA, 
1  EISP,  OMEGA 

COMMON/CC/LB,  FACR 
COMMON/ FF/ A 
COMMON/II/AA,  BB 


P(  -IXt  IX,  0«  lY,  3),  W(  -IXi  IX,  3) 

LB,  FACR 

AC  “IXj  IX^  *“IXs  IX) 

AA(  NDIM  ,  NDIM  ),  BB(  NDIM  ,  NDIM  ) 
B(NDIM) 

C(  -IXi  IX) 

Z(  NDIM  ) 


THETAI,L,X,  PI, 
,T 


SUBROUTINE  SHIP 

THIS  MOVES  THE  PROGRAM  TO  THE  NEXT  TIME  LEVEL 


THE  VALUES  IN  THE  PARAMETER  STATEMENTS  ARE  CONSTANTS 
PARAMETER  (IX  =  20,  lY  =  80) 

PARAMETER  (H  =  .025,  DELTAT  =  .0325) 

PARAMETER  (NDIM  =  2*IX+1) 

C  THE  FOLLOWING  ARE  COMPLEX  ARRAYS 

COMPLEXX16  P(  -IX!  IX,  0!  lY,  3),  W(  -IXi  IX,  3) 
COMPLEX  LB,  FACR 

COMMON/SS/P,  W 

COMMON/TT/KP,KF,INFA,  THETAI,  L,  X,  PI, 

1  EISP,  OMEGA,  T 

COMMON/CC/LB,  FACR 
C 

REAL  KP,  L,  OMEGA,  EISP,  KF,  INFA,  PI,  THETAI 
REAL  X,  T 
INTEGER  I,  J 
C 
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COMMON/FEL/U,V,TP 

tOMMON/KK/IPVT 

COMMON/OO/C 

C 

REAL  KP,  L,  OMEGA,  EISP,  KF,  INFA,  PI,  THETAI 
REAL  X,  T 

INTEGER  IPVT(NDIM),  JOB  ,U,V,  TP 
TP  =  41 
JOB  =  0 

LB  =  CMPLX  (0.0,  1.0) 

LA  =  LXKFXCOS(THETAI) 

RA  =  L»KFKCOS(THETAI) 

C 

DO  42  I  =  -IX+1,  IX-1,  1 

C(I)  =  INFAX(P(I+1,IY,2)+P(I-1,IY,2)+2»P(I,IY-1,2)) 

1  +  (2  -  (4XINFA))«PCI,IY,2) 

42  CONTINUE 
I  =  -IX 

FACR  =  CEXP((2)*(LBXLA)) 

C(I)  =  INFA*(P(I+1,IY,2)+FACR*P(-I-1,IY,2)+2»P(I,IY-1,2)) 
1  +  (2  -  (4XINFA))*P(I,IY,2) 

I  =  IX 

FACR  =  CEXP((-2)X(LBKRA)) 

C(I)  =  INFAX(FACRXP(I-1/IY,2)+P(-I+1,IY,2)+2»P(I,IY-1,2)) 
1  +  (2  -  (4XINFA))*P<I,IY,2) 

C 

DO  10  I  =  1,  NDIM,  1 

B(I)  =  INFAXCCMPLX(0.0,  0.0)) 

10  CONTINUE 
C 
C 

DO  20  I  =  1,NDIM,1 

DO  30  LJ  =  1,  NDIM,  1 

BCD  =  B(I)  +  BB(I,LJ)XP(-1X+LJ-1,IY,1) 

30  CONTINUE 
20  CONTINUE 
C 

DO  150  I  =  1,  NDIM,  1 

B(I)  =  B(I)  +  C(-IX+I-1) 

150  CONTINUE 

CALL  CGESL( AA, NDIM, NDIM, IPVT,B, JOB) 

DO  46  I  =  -IX,  IX,  1 

P(I,IY,3)  =  B(IX+1+I) 

46  CONTINUE 
RETURN 
END 
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DO  19  1  =  -IX,  IX,  1 
DO  22  J  =  0,  lY,  1 
P(I,J,1)  =  P(I,J,2) 

22  CONTINUE 

19  CONTINUE 

DO  20  I  =  -IX,  IX,  1 
DO  23  J  =  0,  lY,  1 
P(I,J,2)  =  P(I,J,3) 

23  CONTINUE 

20  CONTINUE 

DO  29  I  =  -IX,  IX,  1 
W(I,1)  =  W(I,2) 

29  CONTINUE 

DO  <i0  I  =  -IX,  IX,  1 
WCI,2)  =  W<I,3) 

40  CONTINUE 
RETURN 
END 

FUNCTION  FKXP,  TIME) 

PARAMETER  (IX  =  20,  lY  =  80) 

PARAMETER  (H  =  .025,  DELTAT  =  .0325) 

C0MPLEXX16  P(  -IXt  IX,  0i  lY,  3),  W(  -IX.  IX,  3) 
COMPLEXK16  FI,  SI 
COMMON/SS/P,  W 
COMMON/RN/RUNS 

COMMON/TT/KP,KF,INFA,  THETAI,L,  X,  PI, 

1  EISP,  OMEGA,  T 

REAL  KP,  L,  OMEGA,  EISP,  KF,  INFA,  PI,  THETAI 
REAL  X,  T 
REALX8  Q9,Q5,S 


IFCRUNS  .EQ.  1)  THEN 

S  =  0.0 

ELSE 

S  =  TIME+XPXKF»COS(THETAI) 
ENDIF 
09  =  0 

05  =  AMINKTIME,  18.0) 
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Q8  =  DEXP(2<«Q5-6 .0) 
Q7  =  Q8/(Q8+.1/Q8) 
SI  =  CMPLX(Q9,-S) 

FI  =  Q7XCDEXP(S1) 

RETURN 

END 


.  SUBROUTINE  ALPHA 

THIS  SETS  UP  THE  INITIAL  CONDITIONS  OF  THE 
PROBLEM  INCLUDING  THE  DOMAIN  AND  THE 
ARTIFICIAL  BOUNDARY 

THE  VALUES  IN  THE  PARAMETER  STATEMENTS  ARE  CONSTANTS 
PARAMETER  (IX  =  20,  lY  =  80) 

PARAMETER  (H  =  .025,  DELTAT  =  .0325) 

PARAMETER  (NDIM  =  2«IX+1) 

THE  FOLLOWING  ARE  COMPLEX  ARRAYS 
C0MPLEXX16  P<  -IXi  IX,  0i  lY,  3),  W(  -IXi  IX,  3) 

COMPLEX  LB,  FACR 

COMPLEX  A(  -IXiIX,  -IX«  IX) 

COMPLEX  AA(  NDIM  ,  NDIM  ),  BB(  NDIM  ,  NDIM  ) 

COMPLEX  AL(  -40 1  40) 


COMPLEX  LAND 

COMHON/SS/P,  W 

COMMON/TT/KP,KF,INFA,  THETAI,L,  X,  PI, 

1  EISP,  OMEGA,  T 

COMMON/CC/LB,  FACR 
COMMON/FF/A 

COMMON/GG/GAMMA,  BETTA 
COMMON/HH/MINN,  MAXX 
COMMON/II/AA,  BB 
COMMON/KK/IPVT 
COMMON/RN/RUNS 
COMMOH/FEL/U,V,TP 
COMMON/SN/FREQ,  SPD 

REAL  THETAI,  INFA,  KF,  EISP,  OMEGA,  L,  KP 
REAL  XI,  X,  T 

REAL  GAMMA(-40.  40),  BETTA  (-40i40),  PI 

REAL  SPD,  FREQ 

INTEGER  I,  U,  V,  TP,  MINN,  MAXX,  IPVT(NDIM),  RUNS,  Y 

Y  =  IY/2 
WRITE(2,75) 

75  FORMATdX,'  ') 

HEN  =  RUNSXDELTAT 

DO  10  M  =  MINN,  MAXX,1 
AL(M)  =  CMPLX(0.0,  0.0) 

DO  30  I  =  -IX,  IX,  1 
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30 

10 

32 


111 


XI  ~ 

IF(  ABS(I)  .EQ.  IX)  THEM 
DEL  =  .5 
ELSE 

DEL  =  1.0 
ENDIF 

LAND  =  CEXP(-LB*YXH»BETTA(M)  +  (LB»HEN))  „  ,  „ 

AL(M)  =  AL(M)  +  .5XLANDKHXDEL*CEXPC-LBKGAMMA(M)*XI)»P(I»Y,3) 
CONTINUE 
CONTINUE 
WR1TE( 2»  32) 

FORMAT(ix, 'ALPHA  "N"S  FOR  RADIATING  MODES') 

DO  111  M  =  MINN,  MAXX 
WRITE(2,«)CABSCAL<M)) 

CONTINUE 

RETURN 

END 


SUBROUTINE  CLOSE 


C 


C 


C 


THE  VALUES  IN  THE  PARAMETER  STATEMENTS  ARE  CONSTANTS 
PARAMETER  (IX  =  20,  lY  =  80) 

PARAMETER  (H  =  .025,  DELTA!  =  .0325) 

PARAMETER  (XYX=  AO) 

PARAMETER  (TNT=  1000) 

PARAMETER  (NDIM  =  2*IX+1) 

THE  FOLLOWING  ARE  COMPLEX  ARRAYS 
C0MPLEXX16  P(  -IXi  IX,  Ot  lY,  3),  W(  -IXi  IX,  3) 
COMPLEX  LB,  FACR 

THE  VALUES  IN  THE  COMMON  "BLOCKS"  ARE  VARIABLES 
COMMON/SS/P,  W 

COMMON/TT/KP,KF,INFA,  THETAI,L,  X,  PI, 

1  EISP,  OMEGA  ,T 

COMMON/CC/LB,  FACR 
COMMON/HH/MINN,  MAXX 
COMMON/GF/WGGRF , PGGRF 

REAL  KP,  KF,  INFA,  THETAI,  L,  PI,  EISP,  OMEGA 


X,  T 

XX(  0.  41),  WPLOTC  -IXi  IX),  ZZ(OilOOO) 
PPLOSC  -IX:  IX),  SMALL(IOOO) 

II,  RUNS,  SAM,  IR 
PARAMETER  (SAM  =  1000) 

REAL  WGGRF(  SAM  ),  PGGRF(  SAM  ) 

C 

DO  10  II  =  0,  XYX,  1 
XX(II)  =  II 

10  CONTINUE 
C 

DO  11  IR  =  0,  TNT,  1 
ZZ(IR)  =  IR 

11  CONTINUE 
C 

DO  30  I  =  1,1000,1 
SMALL(I)  =  WGGRF(I) 

30  CONTINUE 

DO  35  I  =  -IX,  IX,  1 

PPLOS(I)  =  CDABS(P(I,IY,3)) 

35  CONTINUE 

DO  45  I  =  -IX,  IX,  1 

WPLOT(I)  =  CDABS(W(I,3)) 

45  CONTINUE 
C 


REAL 

REAL 

REAL 

INTEGER 


61 


OOOOO  O  O  <-)  o  o  ooooo  o  o  ooooo< 


CALL  COMPRS 

CALL  PLOTDCZZ, SMALL, 1000,  .TRUE.  , 'LINLIN* , 'CONV  FACTOR^*, 

1  'CONVERGENCE  GRAPHS', 'RUN  NUMBER  +  2000$', 

2  'CONVERGENCE  VALUED') 

CALL  PLOTD(XX,PPLOS,<il,  .TRUE.  , 'LINLIN' , 'FLUID  PRESSURE^', 

1  'PI/2  LENGTH=1.2  FREQ=  3750$ ', 'HORIZONTAL  POSITION^', 

2  'SCALED  PRESSURES ') 

CALL  PLOTDCXX,WPLOT,<H,  .TRUE.  , 'LINLIN' ,' DISP  AT  3750  HZ$', 

1  'SCALED  LENGTH  OF  1 .2$' , 'SCALED  HORIZONTAL  POSITION^', 

2  'SCALED  DISPLACEMENT  OF  PLATE$') 

CALL  DONEPL 

RETURN 
,  END 

SUBROUTINE  FINN 

THIS  SUBROUTINE  MAKES  SURE  WE  HAVE  REACHED  THE  STEADY  STATE 
SOLUTION  OF  THE  PROBLEM,  IE  'EOUALIBRIUM* 

THE  VALUES  IN  THE  PARAMETER  STATEMENTS  ARE  CONSTANTS 

PARAMETER  (IX  =  20,  lY  =  80) 

PARAMETER  (H  =  .025,  DELTAT  =  .0325) 

THE  FOLLOWING  ARE  COMPLEX  ARRAYS 
C0MPLEX«16  P(  -IXi  IX,  Ot  lY,  3),  W(  -IXt  IX,  3) 

COMMON/SS/P,  W 

COMMON/TT/KP,KF,INFA,  THETAI,  L,  X,  PI, 

1  EISP,  OMEGA  ,T 

COMMON/EE/WGTWO,  WGTHE,  WGTOT,  CON 
COMMON/PT/PGTWO,  PGTHE,  PGTOT 
COMMON/RN/RUNS 
COMMON/GF/WGGRF, PGGRF 

REAL  KP,  L,  OMEGA,  EISP,  KF,  INFA,  PI,  THETAI 
REAL  X,  T 

REALX8  WGTWO, WGTHE, WGTOT, CON 
REALX8  PGTWO, PGTHE, PGTOT 
INTEGER  I,  RUNS,  SAM 
PARAMETER  (SAM  =  1000) 

REAL  WGGRF(  SAM  ),  PGGRF(  SAM  ) 

DO  3<i  I  =  -IX+2,  IX-2,  1 


WGTWO  =  W(I,2)KDC0NJG(W(I,2))+  WGTWO 

WuiHt  =  W(I,3)KDCONJG(W(l,5))+  WGTHE 

WGTOT  =  (WGTHE  -  WGTWO)  /  WGTHE 

34  COHTINUE 


DO  44  I  =  -IX  ,  IX  ,  1 

PGTWO  =  P(I,IY,2)KDC0NJG(P(I,IY,2))+  PGTWO 

PGTHE  =  P(I,IY,3)KDC0NJG(P(I,IY,3))+  PGTHE 

PGTOT  =  (PGTHE  -  PGTWO)  /  PGTHE 

44  CONTINUE 


IF  (RUNS.GT.  1999)  THEN 
WGGRF(RUNS-1999)  =  WGTOT 
PGGRF(RUNS-1999)  =  PGTOT 
ENDIF 

WRITECZ  88) 

88  FORMAT(ix,'WGGRF, PGGRF  IN  FINN') 
WRITE(2,»)  WGGRF(RUNS-1999) 
WRITE(2,»)  PGGRF(RUNS-1999) 

RETURN 

END 
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